!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! 
!! Unified layout for iterative solvers
!!
!! \n
!!
!! Iterative solvers have a common structure, but different specific
!! details. This module provides a unified layout.
!<----------------------------------------------------------------------
module mod_iterativesolvers_base

!-----------------------------------------------------------------------

 use mod_utils, only: &
   t_realtime, my_second

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp
 
 use mod_state_vars, only: &
   mod_state_vars_initialized, &
   c_stv

 use mod_output_control, only: &
   mod_output_control_initialized, &
   elapsed_format

 use mod_linsolver_base, only: &
   mod_linsolver_base_initialized, &
   c_linpb

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_iterativesolvers_base_constructor, &
   mod_iterativesolvers_base_destructor,  &
   mod_iterativesolvers_base_initialized, &
   c_itpb

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 !> Linear iterative solver problem
 !!
 !! This type describes a linear system from the viewpoint of an
 !! iterative Krylov method.
 !!
 !! For Krylov methods, the vectors of the Krylov space are considered
 !! homogeneous to the unknown \f$x\f$ (this is consistent with the
 !! fact that they are <em>preconditioned</em> residuals) and are
 !! represented as \c c_stv objects. The vector space operations in
 !! \f$K(A,r^{(0)})\f$ are thus taken directly from \c c_stv. Notice
 !! that a scalar product must be provided, overriding the dummy
 !! implementation of \c c_stv.
 !!
 !! The parameters \c abstol, \c tolerance, \c nmax and \c rmax can be
 !! changed from one call to \c solve to the next one, which is
 !! important, for instance, when building Newton-Krylov solvers.
 !!
 !! \section parallel_execution Parallel execution
 !!
 !! Concerning the function arguments, the input parameters must be
 !! defined on all processors and \c scal must return the correct
 !! value on all the processors, as specified also in \c
 !! mod_state_vars.
 !!
 !! Concerning then the nature of the distributed data arrays, as
 !! discussed in \c c_stv there are two kinds of distributed data:
 !! those for which shared data are assumed to be copies of the same
 !! value and those for which shared data are supposed to be added
 !! together. The solution of the linear system \c x is typically
 !! considered of the first kind while the right-hand side is
 !! typically considered of the second kind. However, since the solver
 !! implementations typically consider the <em>preconditioned</em>
 !! residual, we assume that all the variables \c x, \c r and \c v are
 !! of the first kind, as well as any work array defined by the
 !! solver. It is thus the user's responsibility to ensure that such
 !! quantities are computed consistently in \c pres and \c pkry.
 type, extends(c_linpb), abstract :: c_itpb
  logical :: abstol = .false. !< use absolute tolerance
  real(wp) :: tolerance !< tolerance
  integer :: nmax !< maximum number of iterations
  integer :: rmax !< maximum number of restarts
  procedure(i_itsolver), pointer, nopass :: solver !< solver
  !> MPI rank and communicator
  integer :: mpi_id, mpi_comm
  !> used by the solve subroutine to return the final residual, which
  !! can be useful for the implementation of Newton-Krylov solvers;
  !! <em>notice that this is meaningful only on the root process</em>
  real(wp) :: res
  !> used by the solve subroutine to return the iteration count
  integer :: niter
 contains
  procedure, pass(s) :: factor => itsol_factor
  procedure, pass(s) :: solve  => itsol_solve
  procedure, pass(s) :: clean  => itsol_clean
  procedure, nopass :: working_implementation => it_wi
  !> Preconditioned residual \f$r=P^{-1}(b-Ax)\f$
  procedure(i_pres), deferred, nopass :: pres
  !> Preconditioned Krylov vector \f$v=P^{-1}Ax\f$
  procedure(i_pres), deferred, nopass :: pkry
 end type c_itpb

 !> Iterative solver
 abstract interface
  subroutine i_itsolver(x,sys , maxiterex,niter,res)
   import :: wp, c_itpb, c_stv
   !> linear system
   class(c_itpb), intent(in) :: sys
   !> solution (must contain the initial guess)
   class(c_stv), intent(inout) :: x
   !> true if reached m iterations
   logical, intent(out) :: maxiterex
   !> number of iterations (defined as the number of evaluation of
   !! the lhs which, due to the initial evaluation of the residual,
   !! is equal to the number of Arnoldi iterations +1)
   integer, intent(out) :: niter
   !> residual norms
   real(wp), allocatable, intent(out) :: res(:)
  end subroutine i_itsolver
 end interface

 !> Preconditioned residual
 abstract interface
  subroutine i_pres(r,x)
   import :: c_stv
   implicit none
   class(c_stv), intent(in) :: x
   class(c_stv), intent(inout) :: r
  end subroutine i_pres
 end interface

! Module variables

 logical, protected :: &
   mod_iterativesolvers_base_initialized = .false.

 character(len=*), parameter :: &
   this_mod_name = 'mod_iterativesolvers_base'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_iterativesolvers_base_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
          (mod_kinds_initialized.eqv..false.) .or. &
     (mod_state_vars_initialized.eqv..false.) .or. &
 (mod_output_control_initialized.eqv..false.) .or. &
 (mod_linsolver_base_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_iterativesolvers_base_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_iterativesolvers_base_initialized = .true.
 end subroutine mod_iterativesolvers_base_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_iterativesolvers_base_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_iterativesolvers_base_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_iterativesolvers_base_initialized = .false.
 end subroutine mod_iterativesolvers_base_destructor

!-----------------------------------------------------------------------

  pure function it_wi()
   logical :: it_wi
    it_wi = .true.
  end function it_wi

!-----------------------------------------------------------------------
 
 !> Nothing to do
 !!
 !! In the present implementation, we don't do any preliminary
 !! operation for the iterative solvers. If one is required (such
 !! as factorizing a preconditioner), a status field can be included
 !! in the type which defines the problem and used to make the
 !! required operations whenever required.
 subroutine itsol_factor(s,phase)
  class(c_itpb), intent(inout) :: s
  character(len=*), intent(in), optional :: phase
   ! Nothing to do
 end subroutine itsol_factor

!-----------------------------------------------------------------------

 !> Call the selected solver
 subroutine itsol_solve(x,s)
  class(c_itpb), intent(inout) :: s
  class(c_stv),  intent(inout) :: x
   
  logical, parameter :: verbose = .true.
  logical :: maxiterex
  integer :: niter
  real(wp), allocatable :: res(:)
  real(t_realtime) :: t0, t1
  character(len=50+s%nmax*s%rmax*9) :: message(6)

  character(len=*), parameter :: &
    this_sub_name = 'itsol_solve'

   t0 = my_second()

   call s%solver(x,s , maxiterex,niter,res)

   t1 = my_second()
   if(verbose.and.(s%mpi_id.eq.0)) then
     write(message(1),'(a,i7,a)') '(a,',niter,'e9.1)'
     write(message(2),'(a)') 'Iterative solver diagnostics:'
     write(message(3),'(a,l2)')         '  maxiter:  ',maxiterex
     write(message(4),'(a,i7)')         '  niter:    ',niter
     write(message(5),trim(message(1))) '  residuals:',res
     write(message(6),elapsed_format)   '  el. time: ',t1-t0
     call info(this_mod_name,this_sub_name,message(2:6))
   endif
   if(s%mpi_id.eq.0) then
     s%res = res(size(res))
     s%niter = niter
   else
     s%res = -1.0_wp ! nonsensical value
     s%niter = -1    ! nonsensical value
   endif
   if(allocated(res)) deallocate(res) ! only the root process
 end subroutine itsol_solve

!-----------------------------------------------------------------------

 !> Nothing to do
 subroutine itsol_clean(s)
  class(c_itpb), intent(inout) :: s
   ! Nothing to do
 end subroutine itsol_clean

!-----------------------------------------------------------------------

end module mod_iterativesolvers_base

